*******		REPLICATION FILES
*******		Climate Variability and Irregular Migration to the European Union 
*******		Global Environmental Change 
*******		Fabien Cottier and Idean Salehyan
*******		Replication: Results of the empirical nalyses shown in the main text
*******		This version: April 2021

* note: stata packages required for running models
* - coefplot
* - crossfold
* - plottig


clear all

set more off
set scheme plottig
cd "path/to/replication/directory"

* load data
use "data_quarter.dta", clear 

* load stata risk scripts
qui do scripts/script_rr_spei_082019

* duplicates data
duplicates list cowc year quarter
duplicates tag cowc year quarter, gen(_dupl)
duplicates drop cowc year quarter, force

* set time series
sort cowc year quarter
gen quarter_int = real(regexs(1)) if regexm(quarter,"([0-9]+)")
gen quarterS=(year-2005)*4+quarter_int
order cowc nationality year quarter quarter_int quarterS
tsset cowc quarterS
 
* gen additional variables
encode continent, gen(contN)
recode contN (5=.)


* generate lag quartely migration variables
by cowc: gen nmigrq_exbalk_1tl = nmigrq_exbalk[_n-1]
by cowc: gen nmigrq_exbalk_2tl = nmigrq_exbalk[_n-2]
by cowc: gen nmigrq_exbalk_3tl = nmigrq_exbalk[_n-3]
by cowc: gen nmigrq_exbalk_4tl = nmigrq_exbalk[_n-4]
by cowc: gen nmigrq_exbalk_ln_1tl = nmigrq_exbalk_ln[_n-1]
by cowc: gen nmigrq_exbalk_ln_2tl = nmigrq_exbalk_ln[_n-2]
by cowc: gen nmigrq_exbalk_ln_3tl = nmigrq_exbalk_ln[_n-3]
by cowc: gen nmigrq_exbalk_ln_4tl = nmigrq_exbalk_ln[_n-4]


* gen dummies variables for sample splitting: low agriculture / high agriculture
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* i.year if migr100_exbalk==1, ///
 fe vce(cluster cowc)

sum agriLab if e(sample) & year==2010, d
scalar agriMed=r(p50)
gen agri_H=0 if e(sample) & year==2010
recode agri_H (0=1) if agriLab-agriMed > 0 & year==2010
by cowc: replace agri_H = agri_H[21] if missing(agri_H)


* generate dummies variable for extreme weather (cut-off 10 percentile / 90 percentile)
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* wmean_speiy /// 
 i.year i.quarter_int if migr100_exbalk==1, ///
 fe vce(cluster cowc)

centile (wmean_speiy) if e(sample), centile (10 90)
gen spei_drought = 0 if wmean_speiy!=.
recode spei_drought (0=1) if wmean_speiy <= r(c_1) 
centile (wmean_speiy) if e(sample), centile (10 90)
gen spei_hrain = 0 if wmean_speiy!=.
recode spei_hrain (0=1) if wmean_speiy >=  r(c_2) 


*******		  Macros

global SPEI_Y0 wmean_speiy

global SPEI_SQ_Y0 c.wmean_speiy##c.wmean_speiy
 
global SPEI_Y0Y2 L(0 4 8).wmean_speiy

global SPEI_SQ_Y0Y2 c.wmean_speiy##c.wmean_speiy cL4.wmean_speiy##cL4.wmean_speiy /// 
  cL8.wmean_speiy##cL8.wmean_speiy


*******		 Table 1


* NULL Model // No SPEI term
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 /// 
	i.year i.quarter_int if migr100_exbalk==1, ///
	fe vce(cluster cowc)
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_*  /// 
	i.year i.quarter_int if e(sample), ///
	fe vce(cluster cowc)
est sto tab1_m0
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse null model     : " CVrmse
* list countries included in sample 
est resto tab1_m0
unique cowc if e(sample)
tabulate nationality contN if e(sample)


* Model 1
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 /// 
	i.year i.quarter_int if migr100_exbalk==1, ///
	fe vce(cluster cowc)
est sto tab1_m1
* F test spei parameters & aic
test wmean_speiy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse
* relative change for a -/+ 0.5 SPEI shock
est resto tab1_m1
lincom _b[wmean_speiy]*-0.5, eform
lincom _b[wmean_speiy]*0.5, eform


* Model 2
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_SQ_Y0 ///
 i.year i.quarter_int if migr100_exbalk==1, ///
 fe vce(cluster cowc)
est sto tab1_m2
* F test spei parameters & aic
test wmean_speiy c.wmean_speiy#c.wmean_speiy
estat ic
* F test additional spei parameters
test c.wmean_speiy#c.wmean_speiy
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse


* Model 3
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0Y2 /// 
	i.year i.quarter_int if migr100_exbalk==1, ///
	fe vce(cluster cowc)
est sto tab1_m3
* F test spei parameters & aic
test L0.wmean_speiy L4.wmean_speiy L8.wmean_speiy
estat ic
* F test additional spei parameters
test L4.wmean_speiy L8.wmean_speiy
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse


* Model 4
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_SQ_Y0Y2 ///
	i.year i.quarter_int if migr100_exbalk==1, ///
	fe vce(cluster cowc)
est sto tab1_m4
* F test spei parameters & aic
test L0.wmean_speiy L4.wmean_speiy L8.wmean_speiy ///
 cL0.wmean_speiy#cL0.wmean_speiy cL4.wmean_speiy#cL4.wmean_speiy cL8.wmean_speiy#cL8.wmean_speiy
estat ic 
* F test additional spei parameters
test L4.wmean_speiy L8.wmean_speiy ///
 cL4.wmean_speiy#cL4.wmean_speiy cL8.wmean_speiy#cL8.wmean_speiy
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse




*******		 Table 2

* NULL Model // agrarian countries
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 /// 
	i.year i.quarter_int if migr100_exbalk==1 & agri_H==1, ///
	fe vce(cluster cowc)
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* /// 
	i.year i.quarter_int if e(sample), ///
	fe vce(cluster cowc)
est sto tab2_m0H
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse null model     : " CVrmse
* list of countries
est resto tab2_m0H
tabulate nationality contN if e(sample)

* NULL Model // non-agrarian countries
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 /// 
	i.year i.quarter_int if migr100_exbalk==1 & agri_H==0, ///
	fe vce(cluster cowc)
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* /// 
	i.year i.quarter_int if e(sample), ///
	fe vce(cluster cowc)
est sto tab2_m0L
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse null model     : " CVrmse
global CVrmsenull_LA=round(CVrmse,0.001)
* list of countries
est resto tab2_m0L
tabulate nationality contN if e(sample)


* Model 5 // agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 ///
  i.year i.quarter_int if migr100_exbalk==1 & agri_H==1, ///
  fe vce(cluster cowc)
est sto tab2_m5
* F test spei parameters & aic
test wmean_speiy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse
* relative change for a -/+ 0.5 SPEI shock
est resto tab2_m5
lincom _b[wmean_speiy]*-0.5, eform
lincom _b[wmean_speiy]*0.5, eform

* Model 6 // non-agrarian countries
* low reliance on agriculture
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 ///
  i.year i.quarter_int  if migr100_exbalk==1 & agri_H==0, ///
  fe vce(cluster cowc)
est sto tab2_m6
* F test spei parameters & aic
test wmean_speiy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse
* relative change for a -/+ 0.5 SPEI shock
est resto tab2_m6
lincom _b[wmean_speiy]*-0.5, eform
lincom _b[wmean_speiy]*0.5, eform


* Model 7 // agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_SQ_Y0 ///
  i.year i.quarter_int  if migr100_exbalk==1 & agri_H==1, ///
  fe vce(cluster cowc)
est sto tab2_m7
* F test spei parameters & aic
test wmean_speiy c.wmean_speiy#c.wmean_speiy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse

* Model 8 // non-agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_SQ_Y0 ///
  i.year i.quarter_int  if migr100_exbalk==1 & agri_H==0, ///
  fe vce(cluster cowc)
est sto tab2_m8
* F test spei parameters & aic
test wmean_speiy c.wmean_speiy#c.wmean_speiy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse



*******		 Table 3


* Model 9 // agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* spei_drought spei_hrain ///
	  i.year i.quarter_int if migr100_exbalk==1 & agri_H==1, fe vce(cluster cowc)
est sto tab3_m9 
* F test spei parameters & aic
test spei_drought spei_hrain
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse
* relative change for a severe weather shock
est resto tab3_m9 
lincom _b[spei_drought], eform
lincom _b[spei_hrain], eform

* Model 10 // non-agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* spei_drought spei_hrain  ///
  i.year i.quarter_int if migr100_exbalk==1 & agri_H==0, fe vce(cluster cowc)
est sto tab3_m10
* F test spei parameters & aic
test spei_drought spei_hrain
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse
* relative change for a severe weather shock
est resto tab3_m10
lincom _b[spei_drought], eform
lincom _b[spei_hrain], eform




******************* Additional analyses reported in the main text



****		 correlation irregular migration and asylum applications

preserve

use "data_annual.dta", clear 

* duplicates data
duplicates list cowc year
duplicates tag cowc year, gen(_dupl)
duplicates drop cowc year, force

* set time series
sort cowc year
order cowc nationality year
tsset cowc year

* generate variables
gen applications_ln = log(applications + 1)
by cowc: gen nmigry_exbalk_ln_1tl = nmigry_exbalk_ln[_n-1]

* correlation coefficient
qui xtreg nmigry_exbalk_ln nmigry_exbalk_ln_1tl wmean_speiy /// 
	i.year, ///
	fe vce(cluster cowc)
corr nmigry_exbalk_ln applications_ln if e(sample)


restore



****		 Stationarity of the dependent variable

* Levin-Lin-Chu panel unit-root test
* Schwert 1989 criterion suggests 8 lags

est resto tab1_m1
xtunitroot llc nmigrq_exbalk_ln if e(sample), lags(aic 8) demean



*******		 Additional analyses: Models 5 and 6

**** Seemingly unrelated regression (SUR)

* generate variables
gen nmigrq_exbalk_ln1=0 if nmigrq_exbalk_ln!=.
replace nmigrq_exbalk_ln1=nmigrq_exbalk_ln if agri_H==1
gen nmigrq_exbalk_ln2=0 if nmigrq_exbalk_ln!=.	
replace nmigrq_exbalk_ln2=nmigrq_exbalk_ln if agri_H==0

gen spei1=0 if wmean_speiy !=.
replace spei1=wmean_speiy if agri_H==1
gen spei2=0 if wmean_speiy!=.
replace spei2=wmean_speiy if agri_H==0

gen cowc1=0
replace cowc1=cowc if agri_H==1
gen cowc2=0
replace cowc2=cowc if agri_H==0

gen year1=0
replace year1=year if agri_H==1
gen year2=0
replace year2=year if agri_H==0

gen quart1=0
replace quart1=quarter_int if agri_H==1
gen quart2=0
replace quart2=quarter_int if agri_H==0

gen nmigrq_exbalk_ln1_1tl=0
replace nmigrq_exbalk_ln1_1tl=nmigrq_exbalk_ln_1tl if agri_H==1
gen nmigrq_exbalk_ln2_1tl=0
replace nmigrq_exbalk_ln2_1tl=nmigrq_exbalk_ln_1tl if agri_H==0

gen nmigrq_exbalk_ln1_2tl=0
replace nmigrq_exbalk_ln1_2tl=nmigrq_exbalk_ln_2tl if agri_H==1
gen nmigrq_exbalk_ln2_2tl=0
replace nmigrq_exbalk_ln2_2tl=nmigrq_exbalk_ln_2tl if agri_H==0

gen nmigrq_exbalk_ln1_3tl=0
replace nmigrq_exbalk_ln1_3tl=nmigrq_exbalk_ln_3tl if agri_H==1
gen nmigrq_exbalk_ln2_3tl=0
replace nmigrq_exbalk_ln2_3tl=nmigrq_exbalk_ln_3tl if agri_H==0

gen nmigrq_exbalk_ln1_4tl=0
replace nmigrq_exbalk_ln1_4tl=nmigrq_exbalk_ln_4tl if agri_H==1
gen nmigrq_exbalk_ln2_4tl=0
replace nmigrq_exbalk_ln2_4tl=nmigrq_exbalk_ln_4tl if agri_H==0

* estimate SUR model
qui reg nmigrq_exbalk_ln1 nmigrq_exbalk_ln1_* spei1 i.cowc1 i.year1 i.quart1 if migr100_exbalk==1 & agri_H==1
est sto sureg1
qui reg nmigrq_exbalk_ln2 nmigrq_exbalk_ln2_* spei2 i.cowc2 i.year2 i.quart2 if migr100_exbalk==1 & agri_H==0
est sto sureg2
suest sureg1 sureg2, vce(cluster cowc)

* t test difference coefficient spei
test spei1 = spei2



**** Fixed-effect regression with interaction term for agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* c.wmean_speiy##i.agri_H ///
 i.year i.quarter_int if migr100_exbalk==1, ///
 fe vce(cluster cowc)
est sto inter
